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We use computer simulations to study the glass transition of dense fluids made of polydisperse, 
repulsive spheres. For hard particles, we vary the volume fraction, 99, and use compressible particles 
to explore finite temperatures, T > 0. In the hard sphere limit, our dynamic data show evidence 
of an avoided mode-coupling singularity near i^mct ~ 0.592, they are consistent with a divergence 
of equilibrium relaxation times occuring at ipo ~ 0.635, but they leave open the existence of a 
finite temperature singularity for compressible spheres at volume fraction ip > Lpo- Using direct 
measurements and a new scaling procedure, we estimate the equilibrium equation of state for the 
hard sphere metastable fiuid up to ^po, where pressure remains finite, suggesting that ipo corresponds 
to an ideal glass transition. We use non-equilibrium protocols to explore glassy states above ipo and 
establish the existence of multiple equations of state for the unequilibrated glass of hard spheres, 
all diverging at different densities in the range £ [0.642,0.664]. Classiness thus results in the 
existence of a continuum of densities where jamming transitions can occur. 

PACS numbers: 05.10.-a, 05.20.Jj, 64.70.Pf 
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I. INTRODUCTION 

Making sharp statements about the existence of gen- 
uine phase transitions underlying the observation of non- 
equilibrium glassy states is a favourite game for workers 
dealing with disordered states of matter, one which can 
easily prompt much controversy. Experimentally, a no- 
table exception is the case of spin glasses for which the 
existence of a finite temperature phase transition is not 
controversial because the location of the spin glass 
transition was consistently determined from theoretically 
motivated dynamic scaling relations 0, also reported 
for some frustrated magnets 0]. Despite the ubiquitous 
observation of non-ergodic disordered states across con- 
densed matter physics [1, 0], the existence of genuine 
glassy phases is in fact very rarely established. 

For many-body particle systems Q, such as molecu- 
lar and colloidal glasses, no such scaling predictions arc 
available or experimentally accessible, and the location 
of glass transitions is often deduced from a limited set of 
measurement using uncontrolled extrapolations For 
molecular glasses, fitting the temperature evolution of 
the viscosity of a large number of materials even over 
more than 15 decades cannot qualitatively discriminate 
theories based on the existence of a finite temperature 
singularity from those suggesting a divergence at zero- 
temperature only 0, B Q • Another well-studied instance 
where glassy behaviour is observed is the hard sphere 
system at thermal equilibrium, which becomes highly vis- 
cous when the packing fraction tp increases. Interestingly, 
this idealized model system can be realized experimen- 
tally using colloidal particles [l^l- However, the determi- 
nation of the location of a critical volume fraction where 



the equilibrium relaxation time diverges is plagued by 
uncertainties similar to the ones encountered in thermal 
glasses, since it relies on the extrapolation of a singu- 
larity from a single set of data obtained at increasing 

V PEP. 

In a recent article [ij] , we studied the glassy dynamics 
of a model of soft repulsive particles where the glass tran- 
sition occurs when either the particle volume fraction, ip, 
or the temperature, T, are varied, see Fig. [T] We have 
discovered the existence of activated dynamic scaling in 
the whole (T, ip) plane constraining the actual functional 
form of the divergence, and allowing a precise location of 
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FIG. 1: Sketch of two possible phase diagrams for harmonic 
spheres, (a) In mean-field replica calculations, an ideal glass 
transition occurs at finite temperature for (p above that of 
point G, the ideal glass transition point for hard spheres. The 
pressure of the equilibrium glass diverges at larger density at 
Glass Close Packing (GCP). (b) No glass transition occurs 
before jamming at point J at T = where the equilibrium 
pressure also diverges, with no glass transition occuring at 
finite temperature at larger density. Alternative phase dia- 
grams such as (a) with G = GCP, or (b) with G ^ J are in 
principle possible. 
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the singularity. In particular, we found that in the limit 
T ^ 0, where the particles become infinitely hard, a dy- 
namic singularity seems to occur at a well-defined critical 
packing fraction, ipQ. We called 'point G" the location of 
this special glass point in the phase diagram shown in 

Fig.m 

Our aim in the present work is to study more precisely 
the nature of point G, extending our studies to thermo- 
dynamic observables, most notably equations of state, 
and to explore how point G fits within existing theoret- 
ical frameworks. By simultaneously studying hard and 
compressible particles, we are again able to put severe 
constraints on the nature of the dynamic divergence at 
point G. The picture which is most consistent with our 
data relies on the existence of a hard sphere glass phase 
above ipQ. It remains possible that the ideal glass tran- 
sition reported here is eventually avoided when a much 
larger set of data becomes available, but we show that 
theories where a transition is absent describe the data 
rather poorly. 

We first introduce the theoretical background guiding 
our analysis (Sec.|ll|. We then define the numerical mod- 
els fScc. IIII[) . and successively analyze dynamic (Sec. IIV[) 
and thermodynamic (Sec. |V| behaviour at thermal equi- 
librium, before exploring non-equilibrium glassy states 
(Sec. ED). 



II. IDEAL GLASSES AND JAMMED STATES 

In this section we introduce theoretical ideas, tools and 
predictions needed to analyze the numerical results pre- 
sented below. 



A. Useful analogies for hard and soft particles 

It is useful to draw more precise analogies between the 
glass transition observed in thermal glasses and in hard 
spheres. In thermal glasses, the natural control param- 
eter is the temperature, T, and observables such a equi- 
librium relaxation times, Ta{T), and the energy density, 
e(r), are measured. In the following, we shall measure 
relaxation times by studying the time decay of the self- 
part of the intermediate scattering function, 
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where Vj (t) represents the position of particle j at time t 
in a system composed of TV particles. The brackets rep- 
resent an ensemble average at thermal equilibrium. In 
practice we define Fs{q,Ta) = 1/e for a fixed wavevector 
g, corresponding roughly to the first peak of the static 
structure factor and detecting particle motion typically 
over the interparticlc distance. For isotropic pairwise in- 
teractions defined by the potential V{r), the energy den- 



sity reads 



(2) 



In hard sphere systems, the natural control parameter 
is the volume fraction, 



(3) 



for particles of diameter a and a number density p = 
N/V, where V is the volume of the sample. The equation 
of state of the hard sphere system is then measured by 
defining the reduced pressure, 
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where /cb is the Boltzmann constant, and P the pressure. 
In numerical work, the pressure can be measured from 
the Virial, 



Z = 1- 
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where w{r) = rV (r). For the hard sphere system it can 
conveniently be reexpressed as 



Z((^) = l-H4(^g((7+), 



(6) 



where g(r = a~^) is the pair correlation fimction mea- 
sured at contact. 

It is interesting to notice that the volume fraction if 
plays a role analogous to the energy density e in molec- 
ular glasses, while the reduced pressure is the thermody- 
namic parameter analogous to temperature. Since glassi- 
ness sets in when ip and Z increase in hard spheres, and 
when T and e decrease in molecular glasses, the analogy 
between both system reads: 
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(7) 



Therefore, a finite temperature singularity in thermal 
glasses translates into a singularity at finite pressure in 
hard spheres, while the limits T — *■ and Z — > oo are 
analogous. In both limits, indeed, no particle motion is 
possible. 

At the dynamical level, a reference law for the relax- 
ation time of molecular glasses is the Arrhenius law, 
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where Too and E are two constants with dimensions of 
time and energy, respectively. Deviations from Arrhe- 
nius behaviour, 'super- Arrhenius relaxation', are often 
interpreted as the signature of the non-trivial cooperative 
nature of glassy dynamics. This physical intuition can in 
fact be rigorously established using concepts drawn from 
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linear response theory [1^, [T^, [ij] ■ Angell has introduced 
the notion of 'fragihty' to quantify these deviations [l^ . 
and suggested to represent the experimental data in an 
Arrhenius representation, log Tq, vs. 1/T. Popular func- 
tional forms to account for deviations from Arrhenius 
behaviour are the Vogel-Fulcher-Tamman (VFT) law, 



and for hard spheres, 



Ta{T) = TooCXp 



A 



(9) 



where the exponent S is usually taken to unity, S — I, 
and the Bassler law, 



ra{T) 



B 



(10) 



where the value a = 2 can be obtained from different 
theoretical perspectives [l^, [13, Hi], [l^l • An obvious dif- 
ference between Eqs. © and (fTU]) is the introduction in 
the former of a special temperature Tq where the relax- 
ation time is predicted to diverge, and below which a true 
glass phase should exist. In the latter, equilibrium could 
in principle be maintained down to T = 0, and no glass 
phase exists. 

Using the dictionary in Eq. ^ we now sec that these 
well-known expressions for thermal glasses translate into 
relations between Tq, and Z. Thus the simplest dynamic 
law for hard sphere, analogous to Arrhenius behaviour, 
reads 



,exp[cZ{ip)], 



(11) 



where c is an adimcnsional constant. Just as the Ar- 
rhenius law stems from considering relaxation in a liquid 
arising from a thermally activated local relaxation over a 
fixed energy barrier, Eq. describes a relaxation aris- 
ing from local fluctuations of the volume, which justifies 
fully our analogy between T and 1/Z. This idea is made 
more precise in free volume approaches, see Sec. Ill Cl be- 
low. 

The VFT and Bassler expression respectively become 



Ta{(p) = TooCXp 



and 



,exp[SZ"(^)]. 



(12) 



(13) 



However, since the experimental control parameter is 
the volume fraction while Z is often hard to access 
experimentally, these expressions are usually recast in 
terms of relations between Tq and (p. This can be con- 
fusing since the equivalence between soft and hard par- 
ticles is then easily lost, as it depends on the explicit be- 
haviour of the equilibrium equation of state, Z = Z{(p). 
Fitting formula using Z or if can only take mathemati- 
cally equivalent forms when the pressure is finite at the 
volume fraction where r„ diverges, if it exists. A well- 
known example is the algebraic power laws predicted by 
the mode-coupling theory [2^ for molecular glasses. 



(15) 



[T — Tmct) 
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Alternatively, an instance where the mathematical equiv- 
alence is clearly lost arises if there exists a 'jamming' 
density (p* where Z{ip) diverges. Z (tp* — v?)"^, be- 
cause then an 'Arrhenius' behaviour in Z as in Eq. pT|) 
becomes similar to a 'super- Arrhenius' VFT form when 
the variable ip is used instead. Similarly, strong and frag- 
ile molecular glass-forming materials would look alike if 
the energy density were used instead of the temperature, 
since both types of system could for instance behave as 
Ta exp(c/(e - Go)). 

This example shows that the fragility of hard sphere 
systems should be evaluated by adapting the Angell plot 
to the pressure variable and by plotting logr^ vs. Z, em- 
phasizing possible deviations from the straight line cor- 
responding to the reference law pT|) j24] . Note that in 
a recent experimental work, the Angell plot was adapted 
to soft colloids using (p as an abscissa, instead of the pres- 
sure [l^ . The above considerations suggest that it would 
be very interesting to reanalyze the behaviour of these 
soft colloids along the lines suggested above in order to 
discuss possible changes of 'fragility'. 



B. Ideal glass transition 

The above expressions for the relaxation time in ther- 
mal and hard sphere glasses can be justified from several 
theoretical approaches. We shall describe those theoret- 
ical frameworks which make specific predictions for both 
hard and soft sphere systems, and are therefore relevant 
to the present work. 

A finite temperature/pressure dynamic singularity as 
in Eqs. ^ and and a true glass phase exist in the 
context of random first order transitions [l^, [131 . In the 
fluid phase before the transition, the dynamics is domi- 
nated by the existence of a large number of metastable 
states, allowing the definition of a finite configurational 
entropy. The dynamic transition coincides with the point 
where the configurational entropy of the fluid vanishes, 
and the transition is accompanied by a jump in the spe- 
cific heat (for thermal glasses) or in the compressibility 
(for hard spheres). The glass phase is characterized by 
a vanishing configurational entropy, and the structure of 
phase space is that of a system with one-step level of 
replica symmetry breaking. 

The existence of an 'ideal' glass transition is an exact 
result for a number of many-body interacting models de- 
fined in mean-field geometries, but remains a conjecture 
for three-dimensional realistic particle models. Interest- 
ingly, it is possible to turn this conjecture into a basis for 
actual microscopic, but approximate, calculations of the 
location of the glass transition and the structure of the 
glass phase [1^ [1^. These approximations have been 
applied to fluids of hard and soft particles. For hard 
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spheres, this approach predicts the existence of a diver- 
gence of Tn at a finite pressure, Zg, at a packing fraction 
[23, Ho, HH. Above Lpo the phase which dominates 
the equihbrium measure is a non-ergodic glassy phase 
characterized by replica symmetry breaking. The equi- 
librium equation of state for the glass phase diverges at a 
packing fraction larger than (^o, called 'Glass Close Pack- 
ing', ifGCPi because it corresponds to the densest possi- 
ble glass with an infinite pressure [131 . Replica calcula- 
tions have not been applied to the system of compressible 
spheres studied below, but they would presumably yield 
a phase diagram as sketched in Fig. [T]-(a), with a glass 
transition line Tq(ip > (po) emerging from point G at fo- 
It would be interesting to study also the finite temper- 
ature fate of the jamming transition at GCP from this 
theoretical perspective. 

As mentioned in the introduction, experimental evi- 
dence relics heavily on extrapolations since equilibrium 
can not be achieved very close to the transition. A quite 
favourable experimental finding is the coincidence, with 
reasonable if not decisive accuracy [s^l, of the extrap- 
olated temperatures for the vanishing of the configura- 
tional entropy and for the divergence of the relaxation 
time obtained in several experimental studies of molec- 
ular glass-formers. For hard particles, this coincidence 
seems to hold in numerical studies [s^] , but the extrapo- 
lations are not unambiguous, and arguably not very con- 
vincing. 



C. Free volume ideas 

Free volume arguments are very popular in the litera- 
ture of the molecular glass transition [3J| , although they 
are easier to understand in the context of hard spheres. 
Free volume predictions for hard spheres exist both for 
the equation of state and for the dynamical behaviour. In 
this approach, there exists a maximal packing fraction, 
(p* , where the equilibrium pressure of the hard sphere 
fluid diverges as [35j : 



1 — ipjip*'' 



(16) 



where d is the dimension of space. The critical packing 
fraction is often called 'random close packing' in the 
literature (36j , but we discuss below its meaning in more 
detail. 

To obtain dynamical predictions, one then assumes 
that for if < if* each particle possesses some amount 
of 'free volume', which can then be used to perform lo- 
cal relaxations. In the traditional approach, one gets the 
'Arrhenius' prediction in Eq. (|11[) predicting that Tq di- 
verges at Lp* as 



exp(cZ) ^ exp[cd(^'^/((^* — ip)]. 



(17) 



In this approach, Z and Tq, diverge together at the pack- 
ing fraction ip* , in contrast with the ideal glass transition 
scenario. 



Using Eq. (I7|), we remark that these free volume pre- 
dictions are then the direct analog of the Arrhenius be- 
haviour for thermal glasses, and should serve as a ref- 
erence basis to analyze hard sphere data, the analog of 
'strong' behaviour for thermal glasses. Physically, it in- 
deed makes sense that the divergence in Z should be 
accompanied by a divergence of Tq. since no relaxation is 
possible when all particles touch each other, just as no 
relaxation occurs at T = for molecular glasses. The free 
volume scenario was very recently revisited in Ref. [sj, 
where ip* was given a specific interpretation in terms of 
the jamming transition occuring at point J [ssj in the 
phase diagram of Fig. [l]-(b), predicting that tp* = ipj. 
We discuss the physical meaning of pij further below. 

In recent work, Schweizer and coworkers have devel- 
oped a statistical approach to describe the dynamics of 
hard spheres at large density [s^. They obtain a self- 
consistent equation of motion for the dynamics of a tracer 
particle which they can solve either numerically, or, in 
some limit, analytically. In the latter case, they predict 
in particular a Bassler law as in Eq. (|13p with a = 2. 
They assume further that the pressure diverges at ran- 
dom close packing as in Eq. (jl6p . producing therefore 
the prediction of a non-trivial 'fragile' behaviour for hard 
spheres, namely 



Ta - exp(BZ^) ~ exp[Bd((^*)7((^* - Lp) 



(18) 



Applying this body of ideas to a system of compress- 
ible spheres, we obtain the phase diagram sketched in 
Fig. [l]-(b). Here, the hard sphere system is ergodic up to 
the maximal volume fraction where the equilibrium pres- 
sure diverges and the system jams, and no ideal glass 
transition occurs. At finite temperature, dynamic arrest 
only occurs in the limit where T — > and a true glass 
phase never exists at finite temperature. This phase di- 
agram is consistent with the physical idea that no glass 
transition can occur for thermal glasses at any finite tem- 
perature, because it is always possible to perform 'local' 
relaxations at a finite energy cost. Equivalently for hard 
spheres, one states that as long as some amount of free 
volume is available to each particle, it is possible to relax 
by appropriately displacing a finite set of particles. This 
simple argument, which can be made more formal [ioj . 
provides evidence that the self-diffusion constant for hard 
spheres does not vanish at finite pressure, but it says 
nothing about collective relaxation timescales, which can 
still diverge. Using (with no justification) this argument 
against collective freezing, one would conclude that the 
glass phase only exists along the T = line, as in Fig. [1]- 
(b). Note that self and collective relaxations often yield 
similar results, suggesting that the formal argument in 
Ref. [13] is not very useful for accessible timescales. 

One could easily imagine scenarios intermediate be- 
tween the two sketches in Fig. [TJ For instance, it is 
sometimes stated that no glass transition can occur for 
thermal glasses at any finite temperature, but that hard 
spheres are somewhat different because of the hard con- 
straint imposed by the potential (4l| . In that case, a 
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glass phase would never exist at finite temperature, but 
the glass line at T = could extend down to a low-density 
limit occuring at a finite pressure, so that G ^ J . Note 
however that since the proof presented in Ref. [4l| ap- 
plies to finite size systems, TV < cxd, it does not directly 
contradict the existence of an equilibrium phase transi- 
tion in the thermodynamic limit N ^ oo, as sketched in 
Fig.[I}(a). 

The final alternative is of course the absence of a glass 
of hard spheres before jamming along the T = line, 
together with a glass transition line emerging from ran- 
dom close packing To{(p f*) = 0. This last scenario 
would actually be consistent with free volume models de- 
veloped for both hard and soft particles which does pre- 
dict a VFT divergence at finite temperature for molecu- 
lar glass-formers where the amount of 'free volume' dis- 
appears [i^, but no such glass transition in the case of 
hard spheres where free volume is available up to infinite 
pressure. 



D. Jamming transitions and random close packing 

The glass transition described above locates state 
points where ergodicity is lost. For assemblies of spher- 
ical particles, another transition occurs at large density 
which can be defined in purely geometric terms^without 
invoking thermal equilibrium or ergodicity [ssi . . This 
jamming transition is closely related to the notion of ran- 
dom close packing, very widely discussed experimentally 
in the context of granular materials [iij . 

Granular experiments provide the simplest description 
of a jamming transition for hard particles. Due to grav- 
ity, a disordered assembly of grains will settle in the 
bottom of a container until contact constraints prevent 
any further displacements. The simplest question that 
can be asked is: what is the volume fraction occupied 
by the grains in this jammed state? Turning to fric- 
tionless spheres without gravity, the problem is therefore 
to produce assemblies of particles that can not be com- 
pressed further without allowing overlaps between parti- 
cles [3^ ll^]. These states are therefore infinite 
pressure hard sphere configurations. Note that the re- 
verse is not true, because not all infinite pressure state 
are jammed [47j . 

Since several models for the dynamics of hard spheres 
predict that and Z diverge together at the same den- 
sity, it is very natural to conjecture that glass and jam- 
ming transitions may be two facets of the same transi- 
tion psj . which could therefore be explained in purely 
geometric terms H^. In that case, random close packing 
would correspond to the diverging point in the equilib- 
rium equation of state for hard spheres, and would in 
effect control the glassy dynamics of colloidal particles. 

The situation gets more complicated if an ideal glass 
transition occurs [13, [HO, [5l| because the equilibrium 
pressure does not diverge at the density where ergod- 
icity is lost, and so jamming and glass transitions are 



decoupled, as sketched in Fig. [T]-(a). In that case, the 
energy /volume fraction analogy in Eq. ([7]) suggests that 
finding the glass close packing density is equivalent to 
finding the energy of the ground state in a disordered 
system with a complex energy landscape [soj . This is 
likely a computionally hard problem. This suggests that 
any 'physical' algorithm used to produce T = or Z = 00 
states will end up in configurations that have an energy 
larger than the true ground state (in thermal glasses), or 
that will jam at a density smaller than the glass close 
packed state (in hard spheres), because they will remain 
trapped within the metastable states that are responsible 
for the glassy behaviour at low temperature/large den- 
sity. 

A final complication is the existence for both hard and 
soft particles of a first order transition towards a crys- 
talline structure, which is ignored in the above-mentioned 
approaches. The existence of the crystal phase might 
render the fluid metastable, and crystallization can in- 
deed occur during experiments or simulations for sys- 
tems with low glass-forming ability. It is however easy 
to bypass crystallization in experiments for a large num- 
ber of molecular systems, and for hard spheres using a 
sufficient amount of size polydispersity. If crystalliza- 
tion does not occur (which can be checked by monitoring 
the structure of the fluid), it is then possible to study 
the metastable fluid phase at 'thermal equilibrium', and 
it is the situation usually considered by theoretical ap- 
proaches to the glass transition. In that case, one applies 
concepts from equilibrium thermodynamics to study the 
system in a reduced part of his phase space corresponding 
to metastable disordered states, the crystal region being 
excluded. 

In the same vein, if jamming is studied by driving 
hard sphere conflgurations out of equilibrium using spe- 
cific 'physical' algorithms for compressions, then again 
crystallization is not a critical issue for sufficient polydis- 
persity. However, it should be kept in mind that a large 
number of crystalline and polycrystallinc configurations 
can be built over a broad range of densities encompassing 
the putative glass phase, for instance by artificially mix- 
ing crystalline and fluid phases [s^l . The existence of the 
crystal thus makes the definition of a unique ground state 
for thermal glasses, or the maximum packing fraction for 
hard spheres 'ill-defined' [11], because these concepts are 
not made in reference to thermal equilibrium. Below, 
we shall study jamming transitions specifically designing 
algorithms for which crystallization and demixing are un- 
der controU. 

To make this argument useful to analyze the 
(im)possibility of an equilibrium glass transition (s^ . one 
should additionally show that nucleation of these ordered 
states from the disordered metastable fluid is indeed pos- 
sible at thermal equilibrium. Even then, if nucleation 
barriers happen to be extremely large, then the glass 
transition might indeed be cutoff at extremely large re- 
laxation timescalcs by nucleating the crystal phase rather 
than undergoing an ideal glass transition, but the argu- 
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ment is thus not necessarily useful for experimentally ac- 
cessible timescales. 



III. MODELS FOR HARD AND SOFT 
PARTICLES 



lap has survived after 10^ Monte Carlo steps. For these 
non-equilibrium compressions, we have no insurance that 
Monte Carlo and Molecular Dynamics yield comparable 
results. This is not crucial for our purposes because we 
do not wish to perform dynamic measurements at large 
volume fraction. 



We have explored the equilibrium behaviour and glassy 
states in two numerical models for hard and soft repulsive 
spheres, which we now introduce. 



A. Hard spheres 

We use a standard Monte Carlo algorithm [s^ to 
study numerically a 50:50 binary mixture of hard spheres 
with diameters a and 1.4cr, known to efficiently pre- 
vent crystallization [s^ . We work in a three dimensional 
space with periodic boundary conditions, and mainly use 
N = 1000 particles when studying thermal equilibrium. 
No noticeable finite size effects were found in runs with 
N = 8000 particles performed for selected state points. 
We have detected no sign of crystallization or demixing 
between large and small particles in all our simulations, 
some of which having run for more than 10^*^ Monte Carlo 
steps. 

In an elementary move, a particle is chosen at ran- 
dom and assigned a random displacement drawn within 
a cubic box of linear size O.lcr centered around the ori- 
gin. The move is accepted if the hard sphere constraint 
remains satisfied. One Monte Carlo step corresponds to 
N such attempts. Comparison of Monte Carlo dynam- 
ics with more standard Molecular Dynamics simulations 
of glass-forming liquids have provided evidence that slow 
relaxation in dense fluids is insensitive to the choice of 
a microscopic dynamics (s^ . [Hsj . This is by no means a 
trivial result, as both types of dynamics could in principle 
yield widely different results for the dynamic behaviour, 
especially in those cases where collective particle motions 
are believed to play an important role. 

We have also performed non-equilibrium compressions 
to explore glassy states at large volume fraction. In 
that case we have systematically used two system sizes, 
N = 1000 and N = 8000, and the reported data for com- 
pressions are those for the larger system, although they 
are statistically indistinguishable from those obtained for 
the smaller system. To perform a compression, we use 
the following procedure. We start from an equilibrated 
hard sphere fluid configuration at a given Lp. We then 
perform an instantaneous compression of the simulation 
box, which produces overlaps between particles, which 
are removed using the above Monte Carlo algorithm. As 
soon as all overlaps have disappeared, we perform the 
next compression of the system. We adjust the compres- 
sion rate to maintain the number of overlaps after each 
compression below a constant number, 0.025A^. When 
density gets large, it becomes difficult to remove all over- 
laps and we stop the compression when at least one over- 



B. Harmonic spheres 

We use Molecular Dynamics simulations |53| to study 
a system composed of compressible particles j38l [56j in- 
teracting through a pair- wise potential: V(rij) = e(l — 
Tij/cFij)'^ for Tij < (Jij, V{rij) = otherwise. The inter- 



particle distance is r^. 



r 1 and cr,-. 



((7,+aj)/2, 



where and Oi are the position and diameter of par- 
ticle I, respectively. We shall use the term 'harmonic 
spheres' for this system because the repulsive force be- 
tween spheres is linear in their overlap. We use system 
sizes between 500 and 8000 particles, and report the re- 
sults obtained for TV = 1000, for which no finite size 
effects are detected, within numerical accuracy. 

We prevent crystallisation by using the same 50:50 bi- 
nary mixture of spheres of diameter ratio 1.4 as used in 
the hard sphere case. Up to volume fraction tp = 0.846 
we detect no sign of crystallization at all studied tem- 
peratures; at and above <f = 0.924 there was evidence of 
incipient crystallization at the lowest temperatures p^ . 
However, these crystallization effects occur well away 
from the region of interest around point G in the regime 
(^ = 0.55- 0.75. 

We use e as the energy unit, and ytrl/e as time unit, 
with masses set to unity. All dynamical results are ob- 
tained at thermal equilibrium, which has been carefully 
controlled. When temperature is low and density is large, 
we are not able to thermalize. Crystallization and equi- 
librium issues determine the boundaries of the region of 
investigated state points shown in Fig. [2l 

Finally we have performed exploration of glassy states 
of hard spheres using the compressible sphere system 
as follows. We start from an equilibrated configuration 
of harmonic spheres at a given state point (t/J, T) . We 
then rapidly cool the system at constant density down to 
T — > using Molecular Dynamics. We track pressure and 
energy during compression and obtain valid hard sphere 
configurations at T = when energy vanishes and Z 
remains finite and independent of T at sufficiently low 
temperatures. It is equal to the value of the pressure for 
the corresponding hard sphere system. 



IV. EVIDENCE FOR A DYNAMIC 
SINGULARITY FOR HARD SPHERES 

In Refs. [l^, [5§| we analyzed in detail the volume frac- 
tion dependence of the relaxation time for hard spheres, 
while Ref. [lj| contains a discussion of the interplay be- 
tween density and temperature for harmonic spheres. 
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FIG. 2: Phase diagram for harmonic spheres. Open symbols 
represent the investigated state points. The location of glass 
lines obtained by fitting to known functional form (MCT, 
VFT), and using activated scaling (|22|l are shown, while the 
'super-arrrhenius' line corresponds to volume fractions where 
Tct{<p,T) increases faster than an Arrhenius law when T is 
reduced at constant tp. 



Thus, in this section, wc summarize the main conclu- 
sions drawn in these papers, and complement them by 
further analysis of these dynamic measurements. 



A. Fitting to known functional forms 

For both hard and harmonic spheres, we find that alge- 
braic divergences as predicted by mode-coupling theory 
only hold over a restricted time window of about 2 to 
3 decades in the regime comprising the onset of glassy 
dynamics. This is consistent with findings in other sys- 
tems [2^. First, fitting the relaxation time Taij-p) for 
hard spheres to Eq. (|15p yields a mode-coupling singu- 
larity at (/5MCT ~ 0.592 with a critical exponent 7 w 2.6, 
as reported in Ref. [l3| . This is consistent with previous 
analysis of hard sphere systems [13, HI] ■ 

Second, for harmonic spheres, we fitted the tempera- 
ture evolution of Tq, (1^9, T) at constant Lp between Lp = 0.61 
and if = 0.924 to Eq. ([M]) . We thus obtain a mode- 
coupling transition line, Tmct(v), shown in Fig. [51 
Our data are consistent with a fitted MCT temperature 
which vanishes rapidly as Lp decreases towards (/?mct , al- 
though it is difficult to obtain an accurate determination 
of TivicT very close to i^mct- We find also that the crit- 
ical exponent 7 in Eq. p4[) has a strong volume fraction 
dependence, increasing from 7 w 2.8 for ip = 0.612 to a 
value 7 K, 5.3 for ip ~ 0.8 and above. Although the alge- 
braic divergence predicted by MCT is eventually avoided, 
as we shall describe shortly, it would be very interest- 
ing to analyze the dynamic behaviour of the harmonic 
sphere system using mode-coupling theory: does theory 
reproduce the strong density dependence of the critical 
exponent, and does it predict specific scaling properties 
in the vicinity of the hard sphere point at T = and 



V = <PMCT i 

Since deviations from the algebraic divergence pre- 
dicted by MCT are observed at low enough tempera- 
tures, we repeated our data analysis using standard em- 
pirical approaches. First, we have fitted our data for hard 
spheres using a VFT form. 



Ta {^) ^ exp 



A 



[Vo - pY 



(19) 



by analogy with Eq. ([9]). Imposing the standard value 
5 ~ \, one locates a critical volume fraction at yvFT = 
'Pq{5 = 1) w 0.615. However, using 5 as an additional 
fitting parameter, a slightly better fit is obtained for 
5 « 2.2 from which a larger critical packing fraction, 
ipQ{6 ~ 2.2) w 0.635, is estimated. Remarkably, similar 
conclusions hold for experiments performed on colloidal 
hard spheres [H, [5§| . 

We have then fitted our finite temperature data for 
harmonic spheres to the VFT form in Eq. ([S]), imposing 
an exponent iS = 1 at all densities. We find that such a fit 
can describe our data at all volume fractions rather well 
for ip > 0.636, and report the VFT glass line deduced 
from this fitting procedure in Fig. [21 Again it becomes 
difficult to obtain accurate determination of the VFT 
temperature as ip decreases, but the extrapolation of the 
VFT line is nevertheless in good agreement with the hard 
sphere result </3vft ~ 0.615. 

An interesting outcome of the VFT fitting procedure 
is the density dependence obtained for the fitting param- 
eter A in Eq. For (5 = 1, it is convenient to define 
the parameter D = ATq , which serves as an experimental 
tool to quantify the fragility of supercooled hquids p^ . 
We find that D varies strongly with ip and changes from 
D w 55 for ip = 0.636 down to D w 8 for > 0.80. 
In experimental investigations of the dynamics of super- 
cooled liquids [1^, H^l , D is found to decrease similarly 
from D K 60 for Si02, a strong glass- forming material, 
down to Z? 30 for liquids of intermediate fragilities 
such as glycerol or ZnCl2, and to Z) « 10 for fragile liq- 
uids such as orthoterphenyl (the fit to a VFT form is 
obviously more ambiguous for strong glass-forming ma- 
terials with nearly Arrhenius behaviour). Thus, we find 
that the system of harmonic spheres displays a variation 
in kinetic fragility which encompasses the ran ge o bserved 
in experiments, a central claim made in Ref. |14l |. which 
is confirmed by the present analysis. 

Both MCT and VFT fitting formula contain a diver- 
gence at a critical temperature, corresponding to the 
phase diagram sketched in Fig. [T]-(a). However, to ex- 
plore the possibility of a T = glass line at large density, 
as sketched in Fig. [l]-(b), we have also fitted our data to 
a generalized Bassler form, as in Eq. (|10p . using B and a 
as free fitting parameters for each volume fraction. We 
find that Eq. PH)) also describes our data rather well, 
but we must use an exponent a which increases from 1 
near (p « 0.63 up to a w 3.7 at ip = 0.736, while a = 2 
is traditionnally preferred in supercooled liquids [l^ . 
The range of density where a > 1 delimits the density re- 
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gion where 'super- Arrhenius' behaviour is observed, and 
is indicated with a dashed hne in Fig. [51 

Therefore, we conclude cautiously that fitting our data 
to known functional forms shows that above « 0.63, 
the dynamics at constant slows down with T faster 
than an Arrhenius law, and exhibits dynamics typical of 
fragile glass-forming liquids, with a fragility increasing 
dramatically with ip. This fragility increase is accompa- 
nied by a large variation of the MCT critical exponent 
7, of the VFT fitting parameter D, and the exponent 
a in Eq. (jlOp . but these fitting procedures leave open 
the location of the divergence of Tq,, since both finite 
temperature and singularity-free fitting formula can be 
used to describe our data. Similarly, the activated form 
in Eq. ^9]) is clearly favoured by our hard sphere data, 
possibly with a non-trivial exponent ^ > 1, but the lo- 
cation of (^0 has to be extrapolated from the analysis of 
much smaller volume fractions, and must be discussed 
with caution. 
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B. 



Activated scaling near point G 



The weakness of the above analysis was mentioned in 
the introductory lines of this article. When indepen- 
dently fitting a single data set obtained by changing a 
single control parameter, the range of timescales covered 
by simulations (and experiments!) is usually too small 
to discrimate between very different fitting formula. 

In Ref. [3] we suggested to apply ideas from dynamic 
scaling to the data obtained in the whole (tp, T) plane 
of harmonic spheres to gather more precise information 
on the phase diagram, and in particular the location of 
the dynamic singularity at point G. To the best of our 
knowledge, such an analysis using two control parameters 
has no counterpart in the glass transition literature. 

Our main aim is to determine the location Lp^ of point 
G along the T = hard sphere axis starting from the 
following qualitative considerations about the harmonic 
sphere system. For Lp < Lp^^ the dynamics slows down 
when T decreases, but the relaxation time saturates in 
the limit T ^ to a finite value corresponding to the 
hard sphere fluid. For ip > ipo, however, the relaxation 
time should extrapolate to infinity in this limit, by defi- 
nition of ipQ. These two regimes are obviously delimited 
by (p> = (p>Q^ where the system, like Buridan's ass, 'hes- 
itates' forever between these two regimes. These three 
different situations are all included in the following scal- 
ing form 



Tai(p,T) ~ exp 



A 



T 



(20) 



In this expression, F±{x) are scaling functions applying 
to volume fractions above and below (po, respectively. 
We expect therefore that F-{x — > oo) — * 1 to recover 
the hard sphere fluid limit, Eq. (|19p . when T and 
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FIG. 3: Activated dynamic scaling of relaxation timescales for 
ip £ [0.567, 0.736]. The data for p < p>o and ip > po collapse 
on two distinct branches, as described by Eq. (|20|l . Times 
are rescaled by I/x/T so that the T ^ limit coincides with 
hard spheres thermalized at T = 1. The errorbars describe the 
range of values for which acceptable data collapse is obtained; 
pb = 1.3 is fixed using potential energy considerations. 



Moreover, continuity of Tq at finite T when crossing p ~ 
ipQ implies a common limit for both scaling functions in 
the a; ^ limit: F_{x ^ 0) - i^+(x ^ 0) - x^^'/'^, so 
that Taip = po,T) ~ exp{A/T^^'/^). 

Dynamic scaling was recently observed for athcrmal 
jamming transitions (60l. [6l|. but the nature of the hard 
sphere divergence (algebraic instead of exponential) was 
qualitatively different from Eq. (|20p . and the critical den- 
sity appearing in the scaling formula also had a different 
nature, since these data were not collected at thermal 
equilibrium. 

Using Eq. (^0]) . data at all temperatures T, and for 
volume fractions in the range pi G [0.567,0.736] can be 
collapsed onto the two expected scaling branches and the 
best data collapse is obtained for po = 0.635, S = 2.2 and 
fi ~ 1.3. This is shown in Fig. [31 which reproduces the 
data collapse presented in Ref. [IJ]. Remarkably, while 
the exponent d is not very much constrained by the hard 
sphere data alone, a value below S fa 2 cannot be used 
if data collapse is sought for harmonic spheres at finite 
temperatures. Fixing S near 2 thus allows us to estimate 
ipo is a much more precise manner: 



po = 0.635 ± 0.005, 



(21) 



where the errorbars refer to the range of volume fractions 
for which acceptable data collapse is obtained. 

The value of the exponent fi and the form of the scaling 
variable in Eq. (|20)) were discussed in terms of an effective 
hard sphere radius in Ref. |14| . and the consequences 



p> < Pq. Similarly, F^{x oo) — > oo, for p > pQ. 



on the strong volume fraction dependence of the glass 
fragility explored in some detail. 

While the location of point G is well established by 
the scaling analysis suggested by Eq. (|5D)) . the nature 
of the phase diagram for p > po is not completely de- 
termined since it depends on the specific form of the 
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scaling function F+{x) for large values of its argument. 
If F^{x) diverges at a finite value xq, say F^{x) = 



fj,5/2 /('™^'/2 _ rf.fJ./2\5 



7(^0 

line of the form 



) , then Eq. ([20]) would yield a glass 



(22) 



which is shown as a full line in Fig. [2l However, we 
checked that our data along the branch can also be 
described using a non-diverging form, F^{x) = x'^^^'^{l + 
bx^), which would be consistent with a T = glass line 
above ipQ together with the Bassler form in Eq. PH)) with 
a density dependent effective exponent a. Therefore, we 
conclude that our scaling analysis leaves open the exis- 
tence of a finite temperature singularity above tpo. 

We emphasize that our conclusion that the equilibrium 
relaxation times for hard spheres diverge at the volume 
fraction pTjl relies on a demanding scaling analysis of 
a large set of data in the {f, T) phase diagram of har- 
monic spheres, combined with the analysis of the T = 
hard sphere axis over seven decades of relaxation times. 
However, as is unavoidable in this field, we should not 
exclude that a different dynamic regime can be entered 
when relaxation timescales beyond reach of our numerical 
capabilities are added to the analysis, thereby asymptot- 
ically changing the overall picture presented in this work. 
If so, however, this putative new regime would be exper- 
imentally irrelevant for colloidal particles [l^, [5^ . 



V. EQUILIBRIUM PRESSURE DATA 

The above results support the existence of a non-trivial 
divergence of the relaxation time for hard spheres at a 
critical volume fraction (^o • Quite generally, a divergence 
is expected if the equilibrium reduced pressure Z{ip) also 
diverges, because no particle motion is possible in this 
limit, which is analogous to T = for systems with soft 
potentials, cf. Eq. ([7]). In this section, we ask whether ifQ 
coincides with a divergence in the presssure, and whether 
relaxation times can be related to pressure in a direct 
manner. 



A. Results for hard spheres 

We first describe the equilibrium data obtained from 
direct equilibrium simulations of hard spheres, from 
which pressure is measured through Eq. The results 
are shown as filled circles in Fig. [4] from very low volume 
fractions where Z w 1 up to the largest volume fraction 
for which equilibrium could be reached, Lp — 0.597, where 
Z w 25.3. 

As recalled in Sec. Ill CI free volume arguments pre- 
dict a simple form for the divergence of pressure as 
Z ~ [if* — ip)~^. In Fig. |4] we attempt a descrip- 
tion of our equilibrium pressure data for three values of 
tf*. We impose tp* = i^vft = 0.615, the volume frac- 
tion deduced from a VFT fit to the dynamic data, and 
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FIG. 4: Equilibrium pressure for hard spheres obtained from 
direct Monte Carlo simulations (circles), and extended to 
large tp using the scaling behaviour of the harmonic sphere 
system (triangles). The free volume prediction (|16|l is pre- 
sented for three locations of the divergence: at v'vft = 0.615, 
</5o = 0.635, and the 'best' fit = 0.672, and none of them 
accurately describes the data. Instead, the BMCSL equation 
of state from liquid state theory describes the data very well 
over the entire fluid range and leave the pressure finite at </po, 
where the relaxation time diverges. 



Lp* = ipQ = 0.635, our best estimate for the location of 
point G. These fits are clearly inconsistent with the data, 
as they do not even go through any of the data points. 
This directly implies that the free volume prediction in 
Eq. (|17[) incorrectly represents our data. In the same 
vein, our pressure data are inconsistent with a free vol- 
ume divergence of the equilibrium pressure at (^Qj a-nd we 
do not know how to extrapolate Z{p) to obtain a diverg- 
ing pressure at point G. Thus we conclude that point 
G defined from the study of the equilibrium dynamics, 
does not seem to correspond to point ,/ defined from a 
pressure divergence, in contrast with a recent proposal 
based on a percolation approach [37 1. 

If we insist that Eq. must describe at least the last 



data points obtained at large volume fraction we find that 
the value if* ~ 0.672 represents the 'best' compromise, 
as shown in Fig. [TBI Clearly, however, the shape of the 
pressure is not very well reproduced. Therefore, our re- 
sults are in disagreement with those obtained in Ref. [s^ , 
where the equilibrium equation of state for hard spheres 
was fitted using a free volume expression. We believe that 
the discrepancy stems from the fact that non-equilibrium 
pressure data obtained in fast compressions were incor- 
rectly mixed with equilibrium data and included into the 
free volume fit. 

Finally, we show in Fig. [4| that the so-called Boub- 
lik, Mansoori, Carnahan, Starling and Leland equation 
of state (BMCSL) [H, [63, which is the extension to bi- 
nary mixtures of the Carnahan-Starling equation of state 
for monodisperse hard spheres, describes our equilibrium 
data very accurately over the entire fluid range up to 
(f = 0.597. A similarly good agreement is known to oc- 
cur when using the Carnahan-Starling equation of state 



10 




FIG. 5: 'Angell plot' for hard spheres showing the pressure de- 
pendence of log Ta where the ' Arrhenius' behaviour predicted 
by free volume arguments should appear as a straight line. 
Hard spheres instead display 'fragile' behaviour which is not 
well fitted either by a Bassler law, Eq. (|13|l with a = 2, but is 
consistent with a finite pressure singularity. Excellent fits to 
Eq. Gil) with {5 ^ l,Zo = 30.2) and {S = 2,Zo = 34.4) are 
shown. 



applied to monodisperse systems [64|, up to large vol- 
ume fractions. This excellent agreement is a remarkable 
result, since the BMCSL is a reasonable, but somewhat 
empirically derived, equation of state obtained from in- 
tegral equations using uncontrolled approximations. An 
important feature of the BMCSL equation of state, in the 
present context, is that it is a continuous and derivable 
function up to very large volume fraction, and it predicts 
a singularity in the pressure at an unphysically large vol- 
ume fraction, (p = 1. Therefore, the BMCSL equation of 
state predicts neither a glass transition nor a jamming 
transition, and we cannot use it to extrapolate critical 
values of the volume fraction for the hard sphere system. 

In Fig. [5] we present an 'Angell plot for hard spheres', 
that is, we show the pressure evolution of the logarithm 
of the relaxation time. In this plot, the Arrhenius be- 
haviour Eq. pT|) should appear as a straight line, suggest- 
ing a simple mechanism for the glassy dynamics of hard 
spheres. Clearly, the data in Fig. [5] do not follow such a 
simple law, and hard spheres thus behave in a non-trivial 
manner: they are 'fragile', in the precise sense defined in 
Scc. lIIA.1 A s discussed above, this is in contrast with free 
volume 34 1 and percolation-based [s^l predictions for the 
dynamics of hard spheres. 

We must thus turn to 'fragile' predictions for the be- 
haviour of hard spheres as a function of pressure. We 
first test the prediction by Schweizcr and coworkers (39| 
in Eq. ([18]). We find that a fit of logTc with Z'^ is rather 
In fact, a plot of log Tq, vs. does not linearize 



poor. 

the data, so the quadratic fit in Fig. O is somewhat ar- 
bitrary. Moreover, the fit shown in Fig. O obtained by 
focusing on the data at larger ip yield a microscopic at- 
tempt time Too ~ 8 • 10"'^ (in MC step units), which is 
physically much too small since it is more than 4 orders 
of magnitude smaller than the relaxation time obtained 



in the low density limit, Ta{(p ^ 0) w 10^. We conclude 
therefore that the (asymptotic) expression in Eq. p^ . 
although predicting the correct upward curvature in the 
Angell plot in Fig. [51 is not an accurate representation 
of our data. Thus, the near coincidence between the ex- 
ponent 6 « 2.2 in Eq. and the Bassler expression 
in Eq. (|18p is fortuitous. In fact, using Eq. leaving 
a free to take values different from 2, we find that the 
data are best described, at large density, using a w 6, 
suggestive of a pressure dependence of the data which is 
much stronger than the one predicted by Eq. ([T^ . 

Therefore, we explore a final possibility, suggested by 
the pressure data in Fig. [H of an equilibrium pressure 
which actually stays finite when the relaxation time di- 
verges, as in Eq. p^ . As can be seen in Fig.[5l our data 
are indeed well described by a finite pressure singularity, 
although the hard sphere data themselves are equally well 
fitted using (5 = 1 or 5 = 2, as shown in Fig. \5[ This is 
expected since the volume fraction dependence of Tq, for 
hard spheres was also well fitted by Eq. P^ with 6 = 1 
and 6 = 2, the latter yielding a marginally better fit. To 
conclude with the pressure, we recall that a value close to 
6 = 2 was favoured through the analysis of the dynamic 
data for harmonic spheres. Therefore, our best estimate 
for the critical pressure Zq of the hard sphere system is 
obtained from the fit of the pressure with 6 = 2 shown 
in Fig.[5lin the range Z > 10 and yields: 

ZqW 34.4 ±0.4, (23) 

with errorbars as given by the fitting numerical routine. 



B. Pressure results for harmonic spheres 

Just as investigating the dynamics of harmonic spheres 
in the vicinity of point G allowed for an accurate determi- 
nation of the critical density (po for the divergence of the 
relaxation time for hard spheres, wc can use the harmonic 
spheres to confirm the above finding that the equilibrium 
pressure of hard sphere is finite at ipo i and obtain an in- 
dependent determination of Zq, which does not rely on 
the value of other fitting parameters. 

The temperature evolution of the equilibrium pressure 
obtained in harmonic spheres for densities from tp = 
0.581 up to iy9 = 0.846 is shown in the top panel of Fig. [SI 
For low density, the pressure Z{p, T) smoothly converges 
at low temperature to the pressure of the equilibrium 
hard sphere fluid. Indeed we find that Z{p, T ^ 0) can 
be directly compared to the direct hard sphere simula- 
tions up to <y9 0.6. Of course, it becomes harder to get 
to T — > at higher density because thermalization is too 
hard to achieve within our computer capabilities. 

To analyze these pressure data, we wish to repeat the 
analysis performed for relaxation times in Sec. IIVBI We 
have first attempted to collapse the pressure data assum- 
ing that the equilibrium pressure of the hard sphere fluid 
diverges at point G. Wc thus tried to collapse our data 
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FIG. 6: Top: Temperature dependence of the pressure of har- 
monic spheres at different volume fractions. Bottom: rescal- 
ing the pressure by the factor Z{^p, T —> 0) collapses the data 
for all ifi, in agreement with Eq. (|25|l . For ip > ipo the fluid 
is no longer the thermodynamically stable phase at T = 0; 
these data are shown with open symbols. 

using a critical scaling form near LpQ-. 

(24) 

which obviously has the same interpretation as in 
Eq. ((20|) above. To this end, we imposed the values 
for /i and ipQ obtained in Sec. IIVBI and adjusted the 
function f{if) to obtain the 'best' collapse onto two dis- 
tinct branches for (p above and below i^o ■ This approach 
failed, and we were not able to obtain any collapse us- 
ing this procedure, showing that the scaling properties of 
Tol{^, T) and Z{ip, T) close to point G are very different. 

We then made the opposite hypothesis that pressure 
has a smooth behaviour for the fluid of hard spheres ap- 
proaching point G where it stays finite. This hypothesis 
implies that Z((p, T) is always a smooth function of tem- 
perature. In Fig. [HI we provide evidence supporting this 
hypothesis using the following much simpler scaling as- 
sumption: 

Z(^,T)«Z(^,T^0)5(T), (25) 

where g{T) is a simple function of temperature such that 
g{T ^ 0) = 1. To produce the collapse in Fig. [51 we 



adjusted Z{(p,T 0) for each volume fraction to get 
the best collapse of the data. By definition, this limit 
also corresponds to the equilibrium pressure of the hard 
sphere fluid at this volume fraction. 

We report the results for Z{lp,T 0) in Fig. [H where 
they can directly be compared with the direct measure- 
ments performed in Monte Carlo simulations of hard 
spheres. It is obvious that for ip < 0.6 both data sets 
perfectly overlap, confirming the validity of the scaling 
procedure described in Eq. (|25p in this regime. 

Interestingly, the data collapse in Fig. [H] allows us to 
extrapolate the T — > behaviour of Z{ip, T) for volume 
fractions at which equilibrium hard sphere simulations 
are no longer available. Therefore, the scaling factor 
Z{ip,T 0) shown in Fig. 0] allows us to extend the 
fluid equation of state for hard spheres at larger ip. The 
validity of this extrapolation simply relies on the reason- 
able, but probably only approximately correct, hypothe- 
sis that the temperature dependence of Z(ip, T) does not 
depend on ip over the Hmited interval ip G [0.6,0.635], 
which is precisely the physical content of Eq. (pS)) . 

Moreover, since the equilibrium relaxation timescales 
diverges at (p>o, extrapolation of the fluid branch for ip > 
ipo cannot be performed because a glass phase is present. 
Equilibrium is lost either at a finite T as in Fig. [T]-(a), 
or because the T ^ limit is singular, as in Fig. [T]-(b), 
and so the extrapolated value of the fluid pressure does 
not coincide with the pressure of the thermodynamically 
stable phase above ipo- 

Remarkably, the extension of the fluid branch up to ipo 
shown in Fig. [5] continues to follow the BMCSL equation 
of state. This result was not anticipated, since this equa- 
tion of state is constructed to reproduce the behaviour 
at moderate volume fraction only. As mentioned above, 
the BMCSL pressure only diverges at p = 1, which is 
sometimes interpreted as a weakness of this theoretical 
approach because i^j = 1 is clearly physically accessible. 
Our results suggest the alternative interesting interpre- 
tation that the BMCSL pressure indeed represents the 
equation of state of the fluid up to large volume fractions, 
but the fluid is no longer the thermodynamically stable 
phase above fo — it is a glass [131 . Thus, the divergence 
of the BMCSL fluid equation of state at = 1, and its 
behaviour above ipo are in fact irrelevant. We shall nu- 
merically explore the equation of state of the glass phase 
in the next section. 

Finally, the pressure data obtained using the scaling 
behaviour of harmonic spheres provide an independent 
means to estimate the pressure of the hard sphere system 
at the critical volume fraction ipo. The BMCSL equation 
of state, which represents the scaled data accurately up 
to ipo hits this critical density at the value 

Z{^ = ipo) = 34.4. (26) 

Obviously, the agreement with the independent estimate 
of Zq obtained above in Eq. (|23p is excellent, and provides 
further confidence that the scaling behaviour in Eq. 
can accurately be used. 
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To conclude this section, we have obtained sohd evi- 
dence that the equihbrium pressure of the fluid of hard 
spheres is finite at point G. Thus the phase diagram 
sketched in Fig.[T]-(b) does not hold, and a diverging pres- 
sure occuring when all particles are at contact cannot be 
reached at thermal equilibrium, suggesting that the glass 
transition at point G does not also correspond to a jam- 
ming transition. In other words, wc find that no jam- 
ming transition can be observed at thermal equilibrium. 
In the final section, we shall explore the consequences of 
this finding for understanding the jamming phenomenon. 

VI. EXPLORING MULTIPLE GLASSY STATES 

In this final section, we leave the realm of thermal equi- 
librium to explore glassy states above i^Oj both for hard 
and harmonic spheres. As explained in Sec. IIIII wc shall 
use hard sphere compressions and harmonic sphere an- 
nealings to reach glassy hard sphere states. Our two main 
aims are to investigate the equation of state of glasses 
above Lp^^ and the possibility to observe jamming transi- 
tions with a diverging pressure, which, wc concluded in 
Sec. |Vl cannot be explored at thermal equilibrium. Since 
we abandon thermal equilibrium, we must carefully dis- 
cuss our numerical protocols, because history now be- 
comes part of the story [6^1 • 

A. Choice of non-equilibrium protocols 

To determine the equilibrium equation of state above 
(po, one should in principle get to volume fractions ip > ipo 
while maintaining thermal equilibrium. Since the system 
is non-crgodic, this is not possible in computer simula- 
tions, but can be done in theoretical calculations [28j . 

An intuitive numerical solution could be to compress a 
hard sphere system at a finite compression rate, F, with 
the hope that the limit F ^ can be reached (66l . [67j . 
However, this solution has two immediate drawbacks. 
First, even changing F by a few orders of magnitude, as 
can be done with present day computers, the system falls 
out of equilibrium much above (fio, so that some extrapo- 
lation is again needed [13] . A second problem stems from 
the difficulty in such a non-equilibrium path to check that 
the system is not undergoing some form of crystallization 
(or demixing for a mixture) while being compressed, in 
which case the system could end up in configurations that 
are not necessarily representative of the glass states one 
seeks to investigate. Since we dedicated much effort to 
tackle the ordering issue while studying thermal equilib- 
rium, wc must be similarly careful when studying glasses. 

To circumvent the first of these difficulties, we decided 
to present equations of state obtained during compres- 
sions without attempting any sort of extrapolations. Do- 
ing so, we obtain pressure measurements at large density 
that are upper bound to the true equilibrium pressure, 
since non-equilibrium pressure are larger than the equi- 



librium ones, just as the energy of an annealed glass is 
larger than the equilibrium energy, recall Eq. ([7]). There- 
fore, wc will not be able to investigate the nature of the 
thermodynamic transition at tpo, and the possibility for 
the equilibrium compressibility to have a jump. 

To prevent the exploration of partially ordered or 
dcmixed states, we start our compressions (for hard 
spheres) or annealing (for harmonic spheres) from con- 
figurations that were produced during our exploration of 
the fiuid at thermal equilibrium, for which no tendency 
to order was detected, even in very long simulations. We 
then increase the volume fraction with a very fast com- 
pression rate, or decrease the temperature rapidly, as de- 
scribed in Sec. IIIII We have checked that the particles 
displacements during these very fast compressions arc 
very small, typically much smaller than a particle diam- 
eter, so that our final configurations are no more ordered 
than the original fluid states. Therefore, those states are 
close in spirit to inherent structures usually studied in 
the context of soft potentials [g^. Interestingly, we flnd 
that starting from equilibrium conflgurations, where re- 
laxation was allowed, the pressure measured during com- 
pression or annealing is extremely weakly dependent on 
system size. We have in fact obtained undistinguishable 
resuhs for N = 1000 and N = 8000 particles. This is in 
contrast with infinitely fast annealing from fully random 
configurations [s^, which are very sensitive to system 
size. 



B. Multiple glasses and jamming densities 

Having fixed compression and annealing rates to very 
large values, we are left with a single control parameter 
for exploring glassy states, namely the location of the 
initial equilibrium configuration in the ((/3, T) phase dia- 
gram, which we now vary. 

We first describe the results obtained during hard 
sphere compressions (along the T = axis) starting from 
different initial densities ^pi in the range ipi G [0.35, 0.596]. 
We follow the evolution of the pressure during compres- 
sion using Eq. In Fig. [7] wc present the results ob- 
tained for N — 8000 particles, averaged over 5 indepen- 
dent initial configurations, to obtain a better statistics. 

As soon as the compression starts, the measured pres- 
sure deviates from the equilibrium equation of state, em- 
phasizing that the compression rate is too fast for the 
system to relax to equilibrium, even when volume frac- 
tion is not very large. In the last stages of the compres- 
sion at large volume fraction, the pressure increases very 
rapidly with (f, and appears to be diverging at some fi- 
nal density, (pf (note that Fig. [7] now uses a logarithmic 
scale for the pressure). We stop our compressions when 
Z « 10^, although we could easily continue compressing 
our system up to much larger pressures. 

A central observation in Fig. [7] is that the compression 
curves have a strong dependence on ip\, which survives 
even in the iV ^ oo limit {N = 1000 and 8000 aheady 
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FIG. 7: Various pressure-volume fraction relations obtained 
from simulations of binary hard spheres. Equilibrium data 
and BMCSL line are as in Fig. [l] The non-equilibrium data 
from hard sphere compressions and soft spheres annealing 
demonstrate the existence of multiple branches, all diverging 
at different densities, and result from the glassy behaviour of 
hard spheres. 



yield consistent results), showing that the system can be 
trapped in glassy states that explicitcly depend on the 
preparation history. Since we have carefully controlled 
the protocol to avoid mixing compressions with partial 
ordering, these multiple equations of states result from 
the non-trivial glassy behaviour of the hard sphere sys- 
tem, and not from a competition between randomness 
and order. 

This discussion equally applies to the terminal densi- 
ties Lpi of these compression branches where pressure di- 
verges. These configurations are the hard sphere analogs 
of the inherent structures obtained at T = in systems 
with soft potentials, and l/<pf plays a role similar to the 
potential energy of inherent structures. The fact that 
(y9f has a strong dependence on ipi thus corresponds to 
the observation that the energy of inherent structures 
decreases when the temperature is decreased [6^, as al- 
ready discussed in Refs. [zO, [zll in the context of hard 
spheres. We find that ipf increases slowly from « 0.642 
when ipi is in the dilute regime, but it starts increas- 
ing more rapidly when ipi gets larger than some 'onset' 
volume fraction, (^onsot ~ 0.56, which marks the onset 
of slow dynamics for the hard sphere fluid [7l| . It then 
grows markedly up to (pf « 0.662 for the largest (pi con- 
sidered in this work. Thus our results suggest that glassi- 
ness alone can be responsible for the existence of a finite 
range of volume fractions where hard sphere configura- 
tions get jammed, and give an alternative explanation as 
to why the concept of random close packing is not 'well- 
defined', for which the existence of the crystalline phase 
is irrelevant (isj . 

A final conclusion drawn from the compression curves 
m Fig. [7] is that the pressure at ipQ for non-equilibrium 
compressions is finite. Since these non-equilibrium mea- 
surements represent upper bounds to the equilibrium 
pressure, we obtain a direct verification that the equi- 



librium pressure of the fiuid does not diverge at Lp^ . This 
observation does not involve fits or extrapolations. 

Turning finally to harmonic spheres, we have found 
very similar results. We find that the T limit of the 
pressure obtained in fast annealing of equilibrated config- 
urations at finite temperature also strongly depends on 
the initial temperature, in full agreement with studies of 
inherent structures in systems with soft potentials [6^. 
In Fig. [71 we present the result for the T pressure 
obtained when choosing, for each volume fraction, an ini- 
tial configuration corresponding to the lowest tempera- 
ture for which thermal equilibrium had been reached (see 
the phase diagram in Fig. [2]). The glass equation of state 
obtained this way is qualitatively very similar to the one 
obtained by fast compressions of well equilibrated hard 
spheres, and it diverges near pi « 0.664. 

This rather large volume fraction, (^f — 0.664, pro- 
vides two interesting perspectives. First, if one insists 
that thermal equilibrium can be maintained as long as 
pressure is finite, then one is led to conclude that a glass 
transition with a diverging relaxation time can only oc- 
cur at least above p = 0.664. However, when analyzing 
our equilibrium dynamic data in Sec. IIVI we were never 
able to obtain such a large value for the critical den- 
sity (^0- Second, if we use the theoretical perspective 
sketched in Fig. [T]-(a) with the existence of a true glass 
phase for hard spheres, what we have obtained here is a 
lower bound on the location of the Glass Close Packing 
density, 0.664 < (pccp- 



VII. CONCLUSION 

In this section, we summarize and discuss the main 
results of our numerical study. 

A. Summary 

We studied a binary mixture of hard spheres at in- 
creasing volume fractions. We found that dynamics slows 
down dramatically, and were able to follow the first 7 
decades of relaxation times. We showed that an alge- 
braic power law, as predicted by mode-coupling theory, 
only describes a window of about 3 decades immediately 
after the onset of glassy dynamics near (/Jonsot ~ 0.56. 
The dynamics is best fitted with a generalized VFT law, 
Eq. (|19p . with the best fit obtained for pQ k, 0.635 and 
5 ~ 2.2, but fits with the traditional value 5=1 diverging 
near (/Jvft ~ 0.615 were acceptable. 

We found that the volume fraction dependence of the 
equilibrium pressure was rather modest, and were not 
able to extrapolate these data to obtain a critical vol- 
ume fraction where the pressure diverges, implying that 
Tfj and Z do not appear to diverge at the same volume 
fraction, in contrast with free volume arguments. 

We suggested to build an 'Angell plot for hard spheres', 
representing the evolution of log vs. Z, parametrized 
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by Lp. In this representation, the density analog of the Ar- 
rhenius behaviour, Eq. (jlip . appears as a straight line. 
We find instead that hard spheres display a 'fragile' be- 
haviour, since the pressure dependence of Tq, is much 
more marked. In this representation, we could show that 
several theoretical predictions for the dynamic behaviour 
of hard spheres do not describe our data, the best de- 
scription being offered by the hypothesis that dynamics 
diverges at a finite pressure. We estimated Zq ~ 34.4 for 
the present system. 

We have shown that extending these study using tem- 
perature as a second, independent control parameter, 
could provide an independent confirmation of these con- 
clusions about the dynamic and thermodynamic be- 
haviour of hard spheres. These temperature studies were 
moreover able to make the above statements quantita- 
tively much more convincing. By using dynamic scaling 
ideas, we determined the location (^o of the divergence, 
and the functional form of the relaxation time with a bet- 
ter precision, showing in particular that the most com- 
monly used value (5 = 1 in Eq. (|19p was incompatible with 
our data. For thermodynamic quantities, we discovered 
a much simpler scaling behaviour of the pressure, allow- 
ing to confirm that Z is not divergent at the dynamic 
transition ipi^. 

We have directly confirmed this statement by devising 
non-equilibrium paths to explore glassy states at large 
volume fractions, carefully dealing with crystallization 
and demixing issues. We have found a family of non- 
equilibrium equations of state in the glass phase, which 
can be seen as upper bounds to the true equilibrium 
equation of state above ipQ. We have shown that these 
multiple branches remain distinct in the thermodynamic 
limit; and diverge in a broad range of volume fractions 
between = 0.642 and p = 0.664, where configurations 
get jammed. 



B. Discussion and open questions 

Gathering the above results sheds light on the possible 
phase diagrams for the harmonic sphere system sketched 
in Fig. [U and its T = limit where it coincides with the 
hard sphere system. 

The best description of our data is obtained if we as- 
sume that the relaxation time in the hard sphere limit 
diverges at (/Jq ~ 0.635, where the equilibrium pressure 
of the fluid is « 34.4. Although these numerical values 
are specific to the present binary mixture studied in this 
work, these results invalidate a number of theoretical ap- 
proaches predicting Zq = oo , and more generally the idea 
that the dynamic arrest in hard spheres simply occurs 
when the interparticle distance vanishes [13, [36,, _3X, ■ 
The only theoretical scenario where a glass transition oc- 
curs at finite pressure is the one predicted by theories 
and calculations based on the concept of an ideal glass 
transition of the random first order type [1^, [13, [3l| . So, 
we arc led to the conclusion that our data provide strong 



support for this scenario in the context of hard spheres. 
We were not able to provide a similar evidence for the 
temperature behaviour of harmonic spheres above (poj 
and the possibility of a finite temperature glass line as in 
Fig. [T]-(a) remains an open issue. 

Since we opened the paper with ironic remarks about 
such sharp claims about the existence of glass transitions, 
let us make here a series of cautious remarks about the 
above statement for hard spheres. First, we repeat that 
our conclusions are based upon solid, but necessarily lim- 
ited in range, numerical evidence. Thus, we leave open 
the possiblity that the picture changes when a broader 
range of timescales is covered. But we also insist that the 
range covered numerically is just as large as the range 
covered experimentally in colloids [Tsj , while work in the 
field of molecular glasses suggest that the physics hardly 
changes when several decades of relaxation timescales are 
added. Thus, we can at least claim experimental rele- 
vance for our results. 

Second, although we suggest that the nature of point 
G is consistent with an ideal glass transition of the ran- 
dom first order type, its precise nature remains to be 
established. Theoretically, such a transition is defined 
by the vanishing of the configurational entropy count- 
ing the number of metastable states. Thus, it would 
be important to measure the configurational entropy 
numerically in the vicinity of point G. Unfortunately, 
only approximations to the configurational entropy can 
be accessed numerically, because the very concept of 
metastable states is not well-defined [t^. Another po- 
tential problematic aspect concerns the dynamical be- 
haviour predicted within random first order theory: the 
exponent b « 2.2 in Eq. is usually not the one used 
when analyzing experimental data in molecular glasses, 
and is not the one predicted in Ref. [2^, although scal- 
ing arguments [T^I suggest ways out of the problem that 
remain to be worked out. Also, since explicit replica 
calculations only exist in the hard sphere limit, it would 
certainly be worthwile to extend the computation to har- 
monic spheres at finite temperature. Work is in progress 
in this direction. 

The third remark of caution stems from the existing 
line of research which aims at demonstrating that an 
ideal glass transition cannot exist in hard spheres. In 
Sec. Ill Dl we discussed why we believe published theo- 
retical arguments do not establish the irrelevance of a 
concept of an ideal glass transition in hard spheres, even 
for the specific case considered in Ref. [s^l- Numerical 
evidence against the existence of an ideal glass transition 
was also given in Ref. [T^I- We stress, however, that the 
location of (/jq hi this last work was defined using alge- 
braic laws for the relaxation time, as in Eq. (fT5|) . Using 
a similar description we would have (incorrectly) located 
(po close to v^MCT ~ 0.592, where we indeed were able to 
show that the pressure is not singular, in agreement with 
Ref. [t^I- Therefore, the path to disprove our conclu- 
sions about the nature and location of point G is clear: 
one should provide numerical evidence that thermal equi- 
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librium can be maintained for the present binary mixture 
with no pressure singularity for vohime fractions larger 
than ipo w 0.635 and pressures larger than Zq w 34.4. 
We believe this is a hard numerical task, even when us- 
ing smart algorithms (t^ . 

A final suggestion for future research is the connec- 
tion to jamming transitions suggested in Fig. [71 We 
provided direct evidence that a 'random close packing' 
density can not be provided in a unique manner, as jam- 
ming transitions can occur within a finite range of vol- 
ume fractions. Our argument is different from the one 
presented in Ref. since it does not rely on the ex- 
istence of a crystalline or ordered phase. It is, how- 
ever, in good agreement with recent work using ideas 
from random first order theory to study jamming of hard 
spheres [13, [51} . In future work, we shall study more 
precisely the final configurations obtained when pres- 
sure diverges in Fig. [71 Our conclusion that, along the 
metastable fluid branch, an ideal glass transition inter- 



venes at finite pressure shows that a 'reproducible' jam- 
ming transition cannot be observed at thermal equilib- 
rium. This shows also that the glass transition is likely 
not driven by geometric jamming, but we can tentatively 
use the jamming phase diagram of Ref. [48| in the oppo- 
site direction and claim that glassiness observed in glass- 
forming liquids carries indeed interesting consequences 
for understanding jamming transitions. 
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